Cite this Article
Qin Jia-Xian, Chen Ya-Ming, Deng Xiao-Gang. Stabilized seventh-order dissipative compact scheme for two-dimensional Euler equations. Chinese Physics B, 2019, 28(10): 104701  
Permissions
Stabilized seventh-order dissipative compact scheme for two-dimensional Euler equations
Qin Jia-Xian, Chen Ya-Ming †, Deng Xiao-Gang
College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China

 

† Corresponding author. E-mail: chenym-08@163.com

Abstract

We derive in this paper a time stable seventh-order dissipative compact finite difference scheme with simultaneous approximation terms (SATs) for solving two-dimensional Euler equations. To stabilize the scheme, the choice of penalty coefficients for SATs is studied in detail. It is demonstrated that the derived scheme is quite suitable for multi-block problems with different spacial steps. The implementation of the scheme for the case with curvilinear grids is also discussed. Numerical experiments show that the proposed scheme is stable and achieves the design seventh-order convergence rate.

PACS: 47.11.-j;47.11.Bc;02.60.Cb
Keyword:compact scheme;time stability;simultaneous approximation term;interface treatment
1. Introduction

High-order finite difference methods are well suited for simulations of complex physics as they admit high resolution properties and save large amount of computational resources. Typical examples can be found in fluid dynamics.[1, 2] Although the derivation of high-order finite difference schemes in the interior of the computational domain is quite straightforward, boundary closures often need special investigation[3–5] to avoid accuracy and stability issues. However, it is not an easy task to construct suitable high-order boundary closures to ensure stability.[3]

In our previous work,[6] we showed that a globally seventh-order scheme is not time stable when boundary conditions are imposed directly with the traditional injection method. To rectify this issue, we employed simultaneous approximation terms (SATs)[7] to impose boundary conditions weakly for a dissipative compact finite difference scheme, resulting in a time stable method with global accuracy of the seventh order. The method was demonstrated to be effective for solving one-dimensional (1D) problems, including linear advection equations and Euler equations. The aim of this paper is to extend the algorithm to solve two-dimensional (2D) Euler equations.

To this end, we need to make some modifications to the scheme since SATs involve some subtle issues in the 2D case. As will be seen in section 2, around each corner of the 2D computational domain, there exists a small region where the boundary schemes involve boundary values from two different directions, which is different from the 1D case. Therefore, it necessitates to specially discuss the SATs in these regions to ensure stability.

We consider in section 3 the implementation of the scheme for multi-block problems. In realistic computations the whole computational domain is usually divided into several subdomains with different grid finess to save computational resource considerably (compared with global refinement). For multi-block problems, some interface treatment techniques have to be used.[8–11] Here, the implementation of SATs provides a convenient interface treatment method,[12–14] which has been applied successfully to solve various problems by the summation-by-parts (SBP) community.[15–17]

Since curvilinear grids are often used in practice, we also discuss in section 4 the implementation of the scheme to this kind of problems by focussing on two relevant examples. The SATs are still used in this case to ensure stability. Finally, conclusions are drawn in section 5.

2. Extension to the 2D case

Consider the following two-dimensional Euler equations

where is the conservative variable, and are the fluxes in different directions. Here the total energy is set to be with for air.

2.1. Numerical scheme

To extend the seventh-order scheme considered in Ref.[6] (some details can also be found in Appendix A) to solve Eq. (1), it necessitates to consider its quasi-linear form[18]

where the Jacobian matrices

Here , and is the specific enthalpy. In addition, the Jacobian matrices can be diagonalized as and , where represents the sound speed, and read explicitly as

For a computational domain , we denote the discrete coordinates by with and with , where and are the spacial steps in different directions. Let and be the approximations to and , respectively, and similar notations for other fluxes and matrices. Then the interior semi-discrete scheme for Eq. (1) can be written as

where the expressions for and are presented in Appendix A (see also Ref.[6]). For other points near boundaries, boundary values are involved, which are imposed weakly by using SATs. By introducing the notation
the approximations at the four corner regions read
For other points, the following semi-discrete approximations are used:
Here are penalty coefficients that can be tuned to stabilize the scheme, , , , and stand for the well-posed west, east, south, and north boundary conditions, respectively. In addition, with , where the notation represents the frozen value of at the last time step. Similar notations also apply for . It is observed that the added SATs in Eqs. (7)–(14) vanish when the exact solution is switched into the scheme, meaning that the implementation of SAT method does not contribute to the truncation error and thus does not affect the convergence rate of the original scheme.

Intuitively, one may simply choose the penalty coefficients to be same values as the one-dimensional case, i.e., .[6] To test this choice, we take as an example the vortex model from Ref.[19], whose solution is given by

where . Here the parameters , , , and are used. The computational domain is set to be and the time step is chosen to implement a traditional fourth-order Runge–Kutta scheme. If not stated otherwise, this time scheme will be used for other numerical examples presented in this paper. For simplicity, exact solution values are adopted at the boundaries. One can observe from Fig. 1 that at t=4.2 the density solution near the right corners starts to deteriorate, while at t=6.3 the vortex is almost broken. This phenomenon indicates that some modifications for the scheme should be made to dissolve the non-physical error emerging from the corner regions, where the SATs from two boundaries overlap; see Fig. 2(a).

Fig. 1. Density contours of the vortex obtained by the scheme of Eqs. (6)–(14) with penalty coefficients . Here , the time step is chosen to implement the fourth-order Runge–Kutta scheme. (a) At the moment t=4.2. (b) At the moment t=6.3.
Fig. 2. (a) Illustration of the computational domain. (b) The north–west corner area, where the solution points are affected by SATs from both boundaries.

To show how to improve the scheme, we take as an example the solution point at the west–north corner; see Fig. 2(b). The semi-discrete approximation (9) at this point reads

Here the first penalty term imposes weakly the west boundary condition, whereas the second one imposes the north boundary condition. It is noted that the solution point resides on the west boundary, which means that the second penalty term is not so important here. In this regard, we adjust the coefficient of the first penalty coefficient to weaken the effect of the second one, while other coefficients remain unchanged. The same reasoning applies for other similar points. In this paper, we choose the penalty coefficients to be
for the scheme of Eqs. (6)–(14), although the choice may be optimized further.

2.2. Numerical example

Using the scheme of Eqs. (6)–(14) with the penalty coefficients (20), we recalculate the vortex solution of Eqs. (15)–(18). The new numerical results are presented in Fig. 3, showing good preservation of the vortex. For accuracy test, we introduce the following norm and L 2 norm:

The convergence rate is measured by
where e1 and e2 are the corresponding errors for different numbers of grid cells N1 and N2, respectively. The results shown in Table 1 indicate that the expected seventh-order convergence rate is achieved approximately.

Fig. 3. Density contours obtained by the scheme of Eqs. (6)–(14) with the penalty coefficients (20) on a 41 × 81 grid. Here the time step is chosen to implement the fourth-order Runge–Kutta scheme. (a) At the moment t=4.2. (b) At the moment t=6.3. (c) At the moment t = 10.
Table 1.

Numerical test for the vortex solution of Eqs. (15)–(18) at t = 10. Here the time step is chosen to implement the fourth-order Runge–Kutta scheme.

.
3. Multi-block problems

In this section we intend to show that the developed scheme with SATs is well suited for multi-block problems, which are often used for complex configurations.

3.1. Interface treatment

For simplicity, the computational domain is only divided into two subdomains with an interface situated at see Fig. 4. Denote the conservative variables in the left and right subdomains by with and with , respectively. In addition, assume that are the east boundary values for the left subdomain, and are the west boundary values for the right. Then the problem can be solved by using the scheme of Eqs. (6)–(14) separately for each subdomain.

Fig. 4. Density contours of the vortex of Eqs. (15)–(18) on the grid consists of two subdomains with 41 × 21 and 21 × 21 points. Here the time step is chosen to implement the fourth-order Runge–Kutta scheme. (a) At the moment t = 0. (b) At the moment t = 10.

To verify the scheme, we calculate the vortex convection problem of Eqs. (15)–(18) on the computational domain with the split line x=7.5. The grid ratio is set to be 2:1 in the x direction. Figure 4 depicts that the vortex propagates successfully through the interface. One can also see from Table 2 that the convergence rates in both the two subdomains are well preserved.

Table 2.

Numerical test for the vortex solution of Eqs. (15)–(18) on two different subdomains with the grid ratio 2:1. Here the time step is chosen to implement the fourth-order Runge–Kutta scheme.

.
3.2. Junction point

Next we demonstrate that by exchanging information through penalty terms, the interface technique can naturally handle grid partition with junction points. This time we divide the computational domain into four subdomains with two interfaces situated at and , respectively. Here we still apply directly the scheme of Eqs. (6)–(14) to each subdomain with the boundary values at interfaces setting to be those of corresponding adjacent points; see Fig. 5.

Fig. 5. Illustration of information permutation at the junction point.

Once again we consider the vortex convection problem of Eqs. (15)–(18) to verify the scheme. Here the computational domain is chosen with interfaces situated at x = 10 and y = 0. It is observed from Fig. 6 that the vortex propagates through the interfaces smoothly. In addition, the accuracy test presented in Table 3 shows that the scheme achieves its design order of accuracy approximately.

Fig. 6. Density contours on the grid consists of four subdomains with 41 × 41 (U), 41 × 21 (V), 21 × 41 (W), and 21 × 21 (Z) points. Here the time step is chosen to implement the fourth-order Runge–Kutta scheme. (a) At the moment t = 0. (b) At the moment t = 5. (c) At the moment t = 10. (d) At the moment t = 15.
Table 3.

Numerical test for the vortex solution of Eqs. (15)–(18) at time t = 10 on four different subdomains with the grid ratio set to be the same as that presented in Fig. 6. For brevity only the number of solution points in the subdomain U are presented in the table. Here the time step is chosen to implement the fourth-order Runge–Kutta scheme.

.
4. Discussions on curvilinear grids

In this section, we study how to implement the numerical scheme for problems with curvilinear grids, which are often needed in practice for complex configurations. In this case, we need to consider the problem in the computational coordinates , where the Euler equations (1) read

where the Jacobian , , . Here, we apply the difference scheme presented in Appendix A to calculate numerically the metrics , , , , and the Jacobian J. To ensure stability, it still necessitates to use SATs similarly to the case for Eq. (1) on Cartesian grids; see e.g., Eq. (19). For convenience of the reader, some necessary formulae of the Jacobian matrices of the fluxes and are presented in Appendix B to implement the similar scheme as Eqs. (6)–(14).

4.1. Stationary model

The first example is a model governed by the following Euler equations with a source term

where the source term is determined analytically by so that the problem admits the following stationary solution[20]

The numerical test for the stationary model of Eqs. (23)–(27) is shown in the following Table 4.

Table 4.

Numerical test for the stationary model of Eqs. (23)–(27). Here the time step is chosen to implement the implicit Euler scheme.

.

The curvilinear grid[21] considered here (see Fig. 7(a)) is generated from the standard computational domain by

where the parameters , , and are chosen. The implicit Euler scheme is implemented for time-marching until the residue reaches machine zero (see Fig. 7(b) for the final solution).

Fig. 7. (a) Illustration of a 41 × 41 curvilinear grid. (b) Final density contour obtained on the 41 × 41 grid.
4.2. Nozzle flow

Consider a channel flow[22, 23] governed by the Euler equations (1). The physical domain (Fig. 8(a)) is bounded between x=−0.75 and x=0.75 in the x direction, and is bounded by

in the y direction. The grid is produced analytically by the expressions
with ξ and η being partitioned equally.

Fig. 8. A grid of the channel with 61 × 41 solution points and the corresponding numerical result of the density. (a) The 61 × 41 mesh. (b) Final flow field.

The initial flow field is set to be the free stream with , , , , , where represents the sound speed. The left boundary is set to be subsonic inflow while the right is subsonic outflow. Both the top and the bottom are set to be slip-wall boundary conditions. How to impose the slip-wall boundary condition weakly is explained in Appendix C.

We run the calculation on a grid with 61 × 41 solution points until the residue reaches the machine zero. It can be observed from Fig. 8(b) that the obtained density contour is smooth, showing the effectiveness of the scheme. Similar results can also be obtained on other grids with different finess. As the focus here is just to test whether the scheme is applicable for the case with curvilinear grids and wall boundary conditions, the convergence rate will not be discussed, which depends also on the smoothness of grids and the implementation of boundary conditions.

In addition, to compare the results obtained by the modified scheme and the scheme with the same penalty coefficients as Ref.[6], we show in Fig. 9 the density contours and velocity vectors near the bump at t=0.054. It is observed that while the unmodified scheme starts to break up, the modified one still preserves well, demonstrating again that the modified coefficients are necessary to ensure stability.

Fig. 9. Density solutions and velocity vectors near the bump obtained on the 61 × 41 grid at t=0.054. (a) The unmodified result. (b) The modified result.
5. Conclusion

In this paper, we generalized the globally seventh-order dissipative compact scheme with SATs[6] to the 2D Euler equations. The choice of penalty coefficients for SATs was reconsidered to stabilize the scheme. It was shown that the scheme with SATs is very convenient for dealing with multi-block problems with conformal grids. In addition, the implementation of the scheme for the case with curvilinear grids was also discussed, including the slip-wall boundary condition. Various numerical experiments were performed to verify the proposed scheme. The extension to Navier–Stokes equations may be considered in further work.

Appendix A: Spacial discretizations

Here we revisit briefly the spacial schemes presented in Ref.[6]. For the one-dimensional conservation law

partition the computational domain by using solution points and flux points ( ), where N is the number of solution points, and is the spacial step. Let Uj and be approximations to and , respectively. Notations for other variables and matrices introduced in Section 2 are still adopted. The exact left and right boundary values are given by
Then the derivative of flux can be computed by the following formulae: for (i.e., the interior points),
For near left boundary points we use the following difference schemes:
and the schemes for the right can be obtained by symmetry.

While the values Fj involved in the above difference scheme can be computed directly by , the values are calculated by any approximate Riemann solver (in this paper Roe’s Riemann solver is adopted). Here stands for the left and the right values obtained by using the seventh-order dissipative compact interpolation scheme[1]: For the interior flux points ,

where the dissipation constant is used for the left upwind interpolation, and by setting the interpolation becomes right upwind automatically. Due to the limitation of stencil length, the interpolations at near boundary flux points become

Appendix B: Jacobian matrices of fluxes in Eq. (22)

The Jacobian matrices of and in Eq. (22) are

where J is defined in Eq. (22), , , .

Diagonalizing yields , where , ,

with , , , and .

Similarly, we have , where ,

with , and .

Appendix C: Slip-wall boundary condition

Here we introduce for the two-dimensional Euler equations the slip-wall boundary condition imposed through SATs proposed in Ref.[24]. Consider a Cartesian grid and suppose that we are considering a point situated on the wall boundary. Taking Eq. (11) as an example, i.e.,

with . What we need to do is just to set the normal velocity (here on the Cartesian grid refers to u) as zero, i.e., . Then the wall boundary value is constructed by subtracting the normal velocity from (see Ref.[24] for more information). The implementation of the slip-wall condition for the case with curvilinear grids is similar. Hence the description is omitted here.

Reference
[1] Deng X G Jiang Y Mao M L Liu H Y Li S Tu G H 2015 Comput. Fluids 116 29
[2] Jiang Y Mao M L Deng X G Liu H Y 2014 Comput. Fluids 104 73
[3] Deng X G Chen Y M 2018 J. Comput. Phys. 372 80
[4] Wamg Z D Yang J F Wei Y K et al. 2013 Chin. Phys. Lett. 30 094703
[5] Sun D K Jiang D Xiang N et al. 2013 Chin. Phys. Lett. 30 074702
[6] Qin J X Chen Y M Deng X G 2019 Appl. Math. Mech-Engl. 40 823
[7] Carpenter M H Gottlieb D Abarbanel S 1994 J. Comput. Phys. 111 220
[8] Hesthaven J S 1998 SIAM J. Sci. Comput. 20 62
[9] Kopriva D A Kolias J H 1996 J. Comput. Phys. 125 224
[10] Kozdon J E Wilcox L C 2015 SIAM J. Sci. Comput. 38 A923
[11] Sun Y Z Wang Z J Liu Y 2007 Commun. Comput. Phys. 2 310 https://doi.org/10.2514/6.2006-301
[12] Carpenter M H Nordström J Gottlieb D 1999 J. Comput. Phys. 148 341
[13] Carpenter M H Nordström J Gottlieb D 2010 J. Sci. Comput. 45 118
[14] Lundquist T Malan A Nordström J 2018 J. Comput. Phys. 362 49
[15] Gong J Nordström J 2011 J. Comput. Appl. Math. 236 602
[16] Virta K Mattsson K 2014 J. Sci. Comput. 61 90
[17] Wang S Y Virta K Kreiss G 2016 J. Sci. Comput. 68 1002
[18] Toro E 1999 Riemann Solvers and Numerical Method for Fluid Dynamics—A Practical Introduction Berlin Springer pp. 2–5
[19] Fisher T C Carpenter M H Yamaleev N K Frankel S H 2011 J. Comput. Phys. 230 3727
[20] Wang L Yu M 2018 J. Sci. Comput. 75 253
[21] Chen Y M Deng X G 2019 Comput. Fluids 185 13
[22] Casper J Shu C W Atkins H 1994 AIAA J. 32 1970
[23] Dan X Deng X G Chen Y M Dong Y D Wang G X 2016 Comput. Fluids 129 20
[24] Svärd M Nordström J 2008 J. Comput. Phys. 227 4805